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We study the hysteretic evolution of the random field Ising model (RFIM) at T — when the 
magnetization M is controlled externally and the magnetic field H becomes the output variable. 
The dynamics is a simple modification of the single-spin-flip dynamics used in the //-driven situation 
and consists in flipping successively the spins with the largest local field. This allows to perform 
a detailed comparison between the microscopic trajectories followed by the system with the two 
protocols. Simulations are performed on random graphs with connectivity z — 4 (Bethe lattice) 
and on the 3-D cubic lattice. The same internal energy U (M) is found with the two protocols when 
there is no macroscopic avalanche and it does not depend on whether the microscopic states are 
stable or not. On the Bethe lattice, the energy inside the macroscopic avalanche also coincides with 
the one that is computed analytically with the //-driven algorithm along the unstable branch of the 
hysteresis loop. The output field, defined here as AU /AM, exhibits very large fluctuations with the 
magnetization and is not self- averaging. Relation to the experimental situation is discussed. 

PACS numbers: 75.60.Ej, 75.50.Lk, 81.30.Kf, 81.40.Jj 



I. INTRODUCTION 

The random-field Ising model (RFIM) is one of the 
simplest model to study the combined effects of interac- 
tion and disorder in many-body systems. In particular, 
the response of the RFIM to a slowly varying magnetic 
field at zero temperature pj illustrates the athermal dy- 
namical behavior observed in several experimental sys- 
tems in condensed matter physics such as disordered fer- 
romagnets, superconductors, martensitic materials, etc. 
This response is characterized by avalanches and rate- 
independent hysteresis. Recently, the model has also 
been transposed to the context of finance and human 
behavior @]. 

The aim of the present work is to study the T = 
RFIM in a situation that has not been considered so far, 
when one varies the overall magnetization and not the 
magnetic field (which then becomes a derived quantity 
that we will call the "output" field). More generally, 
we want to describe the behavior of athermal systems 
under control of the extensive variable conjugated to the 
intensive force. This concerns for instance the stress- 
strain curves in shape-memory materials that are usually 
obtained by controlling the deformation of the sample 
and measuring the induced stress 3] . One also uses a 
feedback control that imposes a constant variation of the 
magnetic flux in the case of ferromagnets with a very 
steep magnetization curve 0. 

For a system at equilibrium, it is of course equivalent 
to control the force or the conjugated variable: the sys- 
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tern follows a well-defined curve which corresponds to the 
minimum of the energy or the free-energy. This curve 
may be continuous or discontinuous, as is the case at a 
first-order phase transition. The situation is more com- 
plicated when thermal fluctuations are too small to over- 
come the energy barriers and the system remains far from 
thermodynamic equilibrium on the experimental time 
scale. It then follows a metastable, history-dependent 
path, and there is no reason for observing the same be- 
havior with the two protocols. In fact, there is experi- 
mental evidence that hysteresis loops obtained by vary- 
ing extensive variables display bending-back trajectories 
(with a so-called yield point), and large fluctuations in 
the measured force (or field) |4J, |5j . 

In order to simulate this situation with the T = 
RFIM, one needs to introduce a dynamical rule that 
states how to flip the spins as the magnetization is 
changed. There are of course different ways of locally 
minimizing the energy and the choice for the dynamics 
is not unique, even if one imposes a deterministic rule so 
to get the same result when repeating the simulation. In 
this work, we propose to modify the standard single-spin- 
flip dynamics in a minimal way, so that the new dynamics 
may be considered as the "magnetization-driven" version 
of the dynamics used in the field-driven case 6] . The main 
advantage is that there is a close connection between the 
microscopic trajectories followed by the system with the 
two protocols and the results for the macroscopic quan- 
tities (for instance the internal energy) can be readily 
compared, 

Another and more delicate issue concerns the definition 
of the magnetic field as an output variable. The solution 
that we adopt is again very simple but cannot be con- 
sidered as fully satisfactory. In another recent work 7] , a 
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different approach was proposed, extending the study to 
finite temperatures so to define the field as a Lagrange 
multiplier. Comparison between these two approaches is 
discussed below. Part of our study is performed on a 
Bethe lattice with connectivity z = 4 (or, equivalcntly, 
on random graphs with the same connectivity) . This is to 
benefit from the fact that an almost complete analytical 
description is available in the field-driven casepj. loL fiof. 
Comparing our simulation data with these exact results 
will help in understanding the similarities and differences 
between the two protocols. 

In section [nj we review the model in the usual field- 
driven situation and introduce the modifications in the 
dynamics so to describe the magnetization-driven case. 
The simulation results for the Bethe lattice are discussed 
in section lTTTI and those for the 3-D cubic lattice in section 
IIVI We summarize our main findings and conclude in 
section Ivl 



II. MODEL 

The RFIM with single-spin-flip local relaxation dy- 
namics was specifically introduced for studying the H- 
driven situation. It is thus usually formulated from a 
microscopic Hamiltonian Ti that corresponds to the mag- 
netic enthalpy. For the present study, it is convenient to 
first introduce the internal energy hi : 

u = - Y,s l s ~Y, hA w 

where Si = ±1 are spin variables defined on the sites 
i = 1, . . . , N of a lattice and the first sum extends over 
all nearest-neigbor pairs (the coupling constant is taken 
as the energy unit and set to unity). The random fields 
hi are i.i.d. variables sampled from the Gaussian distri- 
bution pih) — exp(— ft, 2 /2<7 2 )/v / 27TCT with standard devi- 
ation a. The enthalpy Ti is then defined as 

H=U- HM (2) 

where M = Si is the overall magnetization. In the 
following, we consider two types of lattice: a 3-D cubic 
lattice and a Bethe lattice with connectivity z — 4. In the 
first case, numerical simulations are performed on finite 
lattices of size N = L x L x L with periodic boundary 
conditions. In the second case, they are performed on 
random graphs with fixed connectivity z — 4 which pro- 
vide a convenient realization of the Bethe lattice in the 
thermodynamic limit. 



A. //-driven dynamics 

The standard //-driven dynamics consists in locally 
minimizing the enthalpy Ti. As the external field H is 
changed, each spin is aligned with its total local field 



fi + H, where 

j/i 

and the summation is over all the z neighbors j of i. A 
configuration {Si} is then (meta) stable when all the spins 
satisfy the condition. 

St = sign(/, + H) (4) 

One usually starts the metastable evolution with H = 
—oo and all spins Si = — 1. H is then increased until 
the total local field vanishes at a certain site. This first 
occurs for the spin with the largest random field, h" lax . 
This spin is then flipped, which in turn changes the lo- 
cal field at the neigbors and may trigger an avalanche of 
other spin flips. The avalanche stops when a new stable 
configuration is reached. H is then increased again until 
a new spin becomes unstable and the evolution continues 
until all the spins flip up. The upper half of the hystere- 
sis loop is obtained in a similar way by decreasing the 
field from +oo to — oo. Note that the external field H is 
kept constant during an avalanche, which corresponds to 
a complete separation of time scales between the driving 
mechanism and the internal relaxation of the system (the 
dynamics is then referred to as "adiabatic" ) . Because 
the interactions are purely ferromagnetic, the dynam- 
ics has also some remarkable properties: it is abelian[l| 
(the order in which unstable spins are flipped during an 
avalanche is irrelevant for determining the final state) 
and it satisfies return-point memory (i . An important 
feature is the existence of a critical amount of disor- 
der a c below which the hysteresis loops are discontin- 
uous in the thermodynamic limit, the jump in the mag- 
netization corresponding to the occurence of a macro- 
scopic avalanche [3. One has a c « 2.2 for the cubic 
lattice[H and a c = 1.781258... for the Bethe lattice 
with connectivity z = 40. 

Fig. shows an example of an //-driven metastable 
evolution on a random graph with connectivity z = 4. 
For the sake of comparison with the M-driven proto- 
col that is introduced in the next section, we plot the 
internal energy per spin u = hi/N as a function of 
the magnetization per spin m — M/N (both quanti- 
ties being parametrized by the external field H). The 
metastable states {Si}i, {Si}2, {Si}^,... visited by the 
dynamics are represented by triangles while the dashed 
lines in-between indicate the avalanches. Note that the 
total number of states when H is varied from —00 to +00 
depends on the disorder strength and on the particular 
realization of the random fields. Typically, there are only 
a few states when a is small (most avalanches are large) 
whereas the number of states approaches its upper limit 
N when a is large. 

Finally, we want to stress that the energetic barriers 
between the metastable states are strictly defined by the 
dynamics. In fact, the very definition of the metastable 
states (i.e. the stability rule (0J) cannot be separated 
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FIG. 1: (Color on line) Comparison between the if-driven and 
M-driven trajectories in the energy-magnetization plane (u = 
Li/N is the internal energy per spin). Data correspond to a 
single disorder realization with a — 2 on a random graph with 
connectivity z — 4 (N = 10 5 ). The triangles represent the 
states visited by the H-driven dynamics which are separated 
by avalanches (dashed lines). The dots are the states visited 
by the M-driven dynamics. 

from the use of the single-spin- flip dynamics. It has been 
shown recently that a slightly better minimization of the 
enthalpy (obtained by allowing also simultaneous flips 
of nearest-neigbor spins) yields much thinner hystere- 
sis loops while not changing the critical behavior of the 
system [l3j. 



B. M-driven dynamics 

We now define an irreversible dynamics for the case 
where the magnetization of the system is changed ex- 
ternally. There is no external field and the potential 
that has to be minimized (at least partially) is the in- 
ternal energy U. Our goal is to generate a sequence of 
states {Si}i, {Si}2, {Si}^,... when M is increased from 
M = -N to M = +N by elementary steps AM = 2. As 
noted in the introduction, we want this dynamics to be as 
close as possible to the single-spin-flip dynamics used in 
the £f-driven case. For instance, we require that the two 
driving mechanisms become equivalent when the spins 
behave independently and the hysteresis vanishes (either 
because the coupling constant is zero or a — > oo). In this 
limit, one must thus flip, for each value of M, the spin 
with the largest random field, h™^. In the general case, 
we propose to use the simplest "extremal" dynamics: the 
spins are flipped one by one (like in the -ff-driven case) 
and, for each value of M, one chooses the spin that most 
decreases or, at least, less increases the internal energy. 
This is the spin with the largest local field, fl nax , and the 
corresponding change in the energy is AU = —2f™- ax . 
After the spin has been flipped, the local fields /, at the 



neigbors are updated and the same rule is applied until 
all spins are flipped. One obtains a different sequence of 
states when starting from M = +N and decreasing the 
magnetization, which yields an hysteresis loop. It may 
be remarked that this new dynamics bears some similar- 
ity with the "extremal" dynamics used in simple models 
of self-organized criticality (see e.g. Ref. 15]). However, 
in the present case, one never reaches a statistically sta- 
tionary state because each spin in the system flips only 
once and m evolves between —1 and +1. 

By construction, the total number of states visited 
by the dynamics is now N and the crucial feature is 
that this sequence of states contains all the iJ-driven 
metastable states as a subsequence. This is due to the 
abelian property of the ii-driven dynamics and can be 
easily understood by noticing that i) the two dynamics 
start with the same initial state (with all Si = — 1 or 
+1), and ii) the spins that are flipped successively within 
the M-driven dynamics are either those which trigger an 
.ff-driven avalanche or those which are involved in this 
avalanche. In other words, the dynamical rule that has 
been chosen generates a sequence of states that are ob- 
tained by flipping in a certain order the spins involved 
in the ii"-driven avalanches. This is illustrated in Fig. 
which shows the sequences of states obtained with the 
two dynamics in the u — m plane. One can see that the 
M-driven trajectory is a sort of random walk that joins 
the metastable states belonging to the iJ-driven trajec- 
tory. 

A more problematic (but separate) issue concerns the 
definition of the output field H associated to the changes 
in the magnetization. As will be discussed below in more 
detail, one difficulty is that many of the states visited by 
the dynamics are not metastable. This means that is not 
possible to find a field that allows for the condition @ to 
be satisfied for all spins. This is because the local field fa 
at some spins down is larger than the local field at some 
spins up (it is easy to see from Eq. (0J that the condition 
for a microscopic configuration {Si} to be metastable at 
some field H is that /™ m , the minimum value of the 
local field among the spins up, is larger than f™ ax , the 
maximum value of the local field among the spins down). 
In this respect, the present situation is totally different 
from the one considered in Ref. where all the states 
obtained with the M-driven dynamics are stable. An 
additional difficulty is that there is no obvious way to 
define an intensive quantity conjugated to M, playing 
the same role as the Lagrange parameter introduced in 
Ref. • The simple solution that we propose is to define 
the field in such a way that the work needed to go from 
the state at M to the state at M + AM is minimal. The 
field thus identifies with the internal force, 

H{m) = AU/AM = -/™ aa; M ■ (5) 

(In order to facilitate the comparison with the iJ-driven 
dynamics we use the same notation for the external and 
the output field. For the metastable states that are com- 
mon to the two dynamics, the field defined by Eq. (J5J 
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is exactly the external field at which these states become 
marginally stable.) 
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FIG. 2: (Color on line) Comparison between the //-driven 
(dashed line) and M-driven (solid line) trajectories in the 
field-magnetization plane. Data correspond to a single dis- 
order realization with a = 2 on a random graph with con- 
nectivity z — 4 (N = 200). The symbols represent the 
metastable states according to the single-spin-flip dynamics 
(some of them do not belong to the //-driven trajectory, as 
discussed in section III A). 

H(m) is not a monotonously increasing function of the 
magnetization. In fact, as can be seen in Fig. [2 in the 
case of a very small system, it strongly fluctuates with to. 
This is a quite different behavior from the one observed 
for the magnetization in the //-driven case. Note how- 
ever that one can easily deduce the //-driven trajectory 
from the M-driven one (when the magnetization is put 
on the horizontal axis): it is just the envelope function 
that tracks the increasing maxima of //(to). 

A direct consequence of Eq.JSJ) is that the M-driven 
dynamics does not yield any dissipation. Indeed, since 
the work HAM is just equal to the variation of the 
internal energy, the area of any closed loop (that goes 
back to the same microscopic state) is zero. This is to 
be contrasted with the situation in the //-driven case in 
which the work is larger than AU inside the avalanches. 
Although the experimental hysteresis loops obtained in 
M-driven conditions have a much smaller area than the 
//-driven loops0,|3, it is not true that the dissipation is 
zero. We shall come back to this issue in sectionlVlwhere 
we discuss some possible modifications in the definition 
of the field so to avoid this "unphysical" feature. 



C. Fluctuations and self-averaging 

The two preceding algorithms allow to simulate a finite 
system for a given realization {hi} of the random fields. 
Comparison with experiments should be performed by 



considering the limit A" — > oo. It is then desirable that 
the results are self-averaging, i.e. that they do not de- 
pend on a particular realization of the disorder in the 
thermodynamic limit. 

In this respect, the situation is different when one is 
controlling the external field H or the magnetization M. 
In the former case, one can divide a macroscopic sys- 
tem into a large number of macroscopic subsystems that 
are all submitted to the same external field. Then, ac- 
cording to a standard argument 0], away from critical- 
ity, the value of the density of any extensive quantity on 
the whole system (for instance the magnetization M (//)) 
is equal to the average of the (independent) values of 
this quantity over the subsystems. According to the 
central limit theorem, this quantity is distributed with 
a Gaussian probability distribution and (strongly) self- 
aver aging [jj]- On the other hand, in the latter case, 
one cannot decompose a system into subsystems having 
the same magnetization and the standard argument does 
not apply. This implies that i) one must carefully study 
the behavior of the sample-to-sample fluctuations of an 
observable as the system size increases so to conclude 
whether or not this observable is self-averaging, and ii) 
one must be cautious in giving a physical meaning to the 
average over disorder. 

In the following, we analyze the self-averaging char- 
acter of an observable X by performing histograms over 
many disorder realizations for a given size N, so to esti- 
mate the probability distribution Pjy(A). We then study 
the behavior of the variance Vx = {X 2 )m — (^0m as N 
increases (here (.)m denotes the average over disorder 
at constant M, which has to be distinguished from {.)h, 
the average over disorder at constant H). X is (strongly) 
self-averaging if Vx ~ 1/N when A" — * oo. 

III. RESULTS FOR THE z = 4 BETHE LATTICE 

In this section we present the numerical results ob- 
tained by simulating the A/-driven dynamics on random 
graphs with fixed connectivity z = 4. Since small loops 
are rare in these graphs (their typical size is of order 
log AT), the results in the large- N limit are expected to 
converge to the results on a true Bethe lattice, i.e. in 
the deep interior of a Cayley tree. For the sake of com- 
pleteness, we first recall the analytical expressions for the 
average magnetization per spin (m)jj and the average in- 
ternal ene rgy per spin (u)u as a function of the external 
field H HEj: 

(m) H = l-2j2 ( Z )p* n (l-P*) z - n Pn (6) 

71=0 ^ ' 

n=0 \ ' 

x[n{l-p n )+a 2 p{z-2n-H)} (7) 
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where the quantity P* is solution of the equation 



n=0 



(8) 



and the functions p n (H) (n = 0,1..., z) are integrals of 
the Gaussian distribution, 



Pn = 



-J(2n-z)-H 



p(h)dh 



(9) 



(Eq. (7) is obtained by summing the different contribu- 
tions to the internal energy computed in Ref.^j-) In the 
following, we shall also use the expression for the proba- 
bility (per spin) that an avalanche is initiated when the 
field is increased from H to H + dH. It is defined as 
G{H)dH with 



G(H) = ( ~ )P* n [l -P*Y 

n—0 



l p(z-2n-H) (10) 



(this expression corresponds to Eq. (14) in Ref.0 with 
x = l). 

When computing P*{H), it is important to take into 
account the fact that Eq. ((HJ) has three real roots in a 
certain range of H below a c . Then, the magnetization 
curve obtained from Eq. © has an S-shape behavior 
as illustrated in Fig. [3] The correct physical solution that 
gives the lower branch of the hysteresis loop corresponds 
to the smallest root. For a certain value of H, this root 
is not real anymore, P*{H) jumps to the largest root 
and there is a discontinuity in the magnetization curve 
associated to the occurence of an infinite avalanche. In 
this context, the intermediate, unstable branch of the S- 
shape curve (from points A to B in the figure) has no 
physical meaning (on the other hand, the branch BC can 
be reached via first-order reversal curves obtained from 
the descending branch of the hysteresis loop, as noted in 
Refs.[ll[l3). 



A. Fraction of stable states along the //-driven and 
A/-driven trajectories 

As can be seen in Fig- El which corresponds to the sim- 
ulation of a small system, there are a few states along the 
M-driven trajectory that are metastable although they 
do not belong to the //-driven subsequence (this could be 
seen as well in the u — m diagram since this property does 
not depend on the definition of the field). Of course, the 
fact that the state visited by the dynamics is stable or not 
depends on the disorder realization. Is is therefore useful 
to introduce the quantity Q{m) that represents the aver- 
age fraction of states that are stable. As shown in Fig. 
these additional metastable states appear less and less 
frequently as the system size increases and there is strong 
numerical evidence that they completely disappear in the 
thermodynamic limit. Therefore, when TV — > oo, Q(m) 
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FIG. 3: (Color on line) Ascending branch of the //-driven 
hysteresis loop on a Bethe lattice with connectivity z = 4 
for g = 1.6. The solid line corresponds to the solution of 
Eqs. JSJ and|HJ. The symbols represent all the metastable 
states obtained along the M-driven trajectories for 10 disorder 
realizations on random graphs of size N = 10 4 (a) and 10 5 
(b). 



also represents the average fraction of metastable states 
along the //-driven trajectory. 

The results of the simulations with the M-driven dy- 
namics below and above a c and different system sizes are 
shown in Fig. For a > <r c , it is found that Q(m) is 
minimum in the range of m that corresponds to the steep- 
est part of the //-driven magnetization curve where the 
avalanches are the largest. For a < <r c , the interval where 
Q(m) = exactly corresponds to the range of the infinite 
avalanche (including in the portion BC of the magnetiza- 
tion curve, as shown in the inset of Fig.0J. The fact that 
Q(m) is strictly smaller than 1 (except for m = ±1) de- 
serves some explanation. With the //-driven algorithm, 
there is indeed a certain probability, for a finite system, 
that a given value of M corresponds either to an hori- 
zontal portion of the magnetization curve (a metastable 
state) or to a vertical jump (an avalanche). Although 
the magnetization curve is continuous in the thermody- 
namic limit (except for the jump below a c ), the probabil- 
ity of "hitting" a metastable states remains smaller than 
1 when N — > oo. In other words, Q(m) tracks the ran- 
dom presence of "holes" in the magnetization curve that 
correspond to the avalanches. This suggests that Q(m) is 
related to the probability of having an avalanche between 
H and H + dH (where H is the field corresponding to m 
in the thermodynamic limit). This probability is given 
by the quantity G(H) defined by Eq. (fTTTf) and the seeked 
relation is 



Q{m) = 2G{H) 



dH 
dm 



(11) 
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slope of the magnetization curve in the thermodynamic 
limit is given by 



FIG. 4: (Color on line) Fraction Q(m) of stable states along 
the ascending branch of the hysteresis loop for a Bethe lattice 
with 2 = 4. The points A, B and C are the same as in Fig. [3 
The symbols are the results of the numerical simulation of 
the M-driven algorithm. The solid line corresponds to the 
analytical expression [11 1 1 and the dashed line indicates the 
infinite avalanche below er c . The inset in (a) shows that Q(m) 
converges to between the points B and C in the thermody- 
namic limit. Simulation data, indicated by different symbols 
are joined by guides to the eye. 



where dH/dm is the inverse slope of the magnetization 
curve in the thermodynamic limit, a quantity that is eas- 
ily computed from Eq. ©. One can see in Fig. 01 that 
the agreement between the simulations and the analytical 
formula is indeed very good. The proof of Eq. re- 
lies on the assumption that the occurrence of avalanches 
(as H is monotonously increasing) corresponds to a non- 
stationary Poisson process 0. For a finite system of size 
N, an avalanche then occurs in the interval dH with a 
rate dP/dH = NG(H) (recall that G{H)dH, as calcu- 
lated in Ref.jjJ, is a probability per spin). The mean 
range of stability (AH)h of a metastable state before an 
avalanche occurs is given by the inverse of the rate, i.e. 



1 



NG(H) 



(12) 



(Note that this quantity becomes infinitesimal in the 
thermodynamic limit.) Since only metastable states con- 
tribute to the variation of H in the interval M, M+2 (the 
field is kept constant during an avalanche), the average 



f = Q(m) hm 

dm y ' n->oo AM 



(13) 



which yields Eq. i|Tl"|) . 

It is interesting to remark that Eq. gives a fi- 

nite, positive value of Q(m) between points B and C in 
Fig- EI a) when one computes m(H) using the largest root 
of Eq. ©. Indeed, as already noted, this part of the H- 
driven hysteresis loop can be reached by an appropriate 
field history starting from saturation: this implies that 
the fraction of metastable states is not zero. On the other 
hand, one gets a meaningless negative value for Q(m) 
between points A and B and the correct physical result 
Q(m) = 0, that states that all configurations are unsta- 
ble, is recovered by setting dm/dH — > oo in Eq. Ulljl . 



B. Internal energy 

We now discuss the simulation results for the internal 
energy per spin u{m). In Fig. [SI we plot on a log-log scale 
the variance V u (m) = (m(to) 2 )m — { u { m )) 2 M f° r selected 
values of to as a function of the system size N. Both 
above and below cr c , it is found that V u (m) decrease like 
1/iV, showing that the energy is a strongly self-averaging 
quantity, a result that is not a priori obvious. Accord- 
ingly, we shall now use the average value of tz(to) to com- 
pare with the exact results obtained with the ii-driven 
dynamics in the thermodynamic limit, although (u(to))m 
may not be a well-defined physical quantity, as remarked 
above. Rather, this must be considered as a convenient 
way of suppressing sample-to-sample fluctuations. 

The comparison is performed in Fig.|S|where (u(ui))h 
is obtained by plotting (u)h as a function of (m)jj, the 
field H being considered as a parameter (see also Fig.^l. 
When a < er c , there is a jump in the magnetization 
and the corresponding discontinuity in (u(to))# is rep- 
resented by a dashed line, the solid line representing the 
internal energy along the intermediate, unstable part of 
the magnetization curve. It can be seen hat the behavior 
of u(m) changes with a. For large disorder, the energy 
has a "double well" structure whereas there is a single 
well when the disorder is small. The change in the be- 
havior occurs at a ~ 2.0 and is therefore not related to 
the critical value of the disorder at which the disconti- 
nuity in the iJ-driven magnetization curve disappears. 
Note moreover that the curves are not symmetric with 
respect to m = 0: there is indeed hysteresis when the 
magnetization is increased from —1 or decreased from 
+1. 

The most remarkable feature in Fig.^is that the av- 
erage internal energy obtained with the M-driven algo- 
rithm appears to coincide with the analytical curve ob- 
tained from Eq. Q in the thermodynamic limit even 
when to is in the range of the macroscopic avalanche for 
(7 < a c (the agreement is better than 10 -3 for N = 10 4 ). 
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FIG. 5: (Color on line) Variance Vu(m) of the internal energy 
per spin it(m) on random graphs with connectivity z = 4 for 
selected values of m as a function of system size. The number 
of disorder realizations is 10 4 . For the sake of clarity, the 
variances for m = and m = 0.5 are divided by 10 and 100, 
respectively. The lines are fits to the form V u (m) ~ N~ p , 
yielding p ~ 1.0 in all cases. 
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FIG. 6: (Color on line) Average internal energy per spin 
on the Bethe lattice with z — 4 below and above a c . The 
symbols are the results of the simulation of the M-driven 
algorithm on random graphs of different sizes with an average 
over 10 4 disorder realizations. The solid line corresponds to 
the analytical expression given by Eq. (7). The dashed line 
indicates the discontinuity associated to the infinite avalanche 
for a < g c . 



Since this is a surprising result, we have carefully checked 
the behavior as a function of the system size (note inci- 
dentally that finite-size effets are not negligible in the 
.ff-driven case as well: this is an issue that has not yet 
been investigated, as far as we know). 

When all avalanches are of microscopic size, the co- 
incidence of the energy along the two trajectories is 
due to the fact that the stable states before and after 
the avalanche (and therefore all the unstable states in- 
between) differ only by a finite (i.e., non-extensive) num- 
ber of spin-flips. Accordingly, the energy of these states 
cannot differ by an extensive quantity and one has 

(u(m)) M = {u{m))t? We = {u{m))Z staUe (14) 

in the thermodynamic limit, as can be checked numeri- 
cally. Moreover, since both the energy and the magneti- 
zation are self-averaging quantities, the averages at fixed 



to or fixed H yield the same result. Therefore, 

(u(m)) H = (u(m))T ble = Hm))T le = (u(m))M ■ 

(15) 

It is more surprising that the equality (it(m))ij = 
(u(m))M is also satisfied inside the macroscopic 
avalanche if one uses the "unphysical" root of Eq. (jHJ 
to compute (u(m))H along the unstable branch of the 
ff-driven magnetization curve. We have no obvious ex- 
planation for this result but we want to stress that it cru- 
cially depends on the order in which the spins are flipped 
during an avalanche. Indeed, as illustrated in Fig. [71 a 
different curve (u(to))m is obtained if one decides for 
instance to flip the spin that less decreases the energy. 
Therefore, the M-driven dynamics that has been chosen 
is precisely the one that yields a gree ment with the ana- 
lytical solution computed in Ref. [l(J- This suggests that 
behind the probabilistic computation in Ref. |£| there is 
perhaps some hidden minimization principle that fixes 
unambiguously the trajectory along the unstable branch. 
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avalanche, one flips the spin that less decreases the energy. 

FIG. 8: (Color on line) Ascending branch of the M-driven 
trajectory in the field-magnetization plane. Data correspond 
to a single disorder realization on a random graph with con- 
C. Statistical behavior of the output field nectivity z = 4 (N — 10 4 ). 



As illustrated by Figs. [21 and [SJ the output field H(m) 
defined by Eq. (5) displays a sporadic, discontinuous be- 
havior with magnetization. When a spin flips, the local 
field at the neigbors is changed by ±2, and this is in- 
deed the approximate size of the fluctuations observed 
in Fig. [SJ This of course does not depend on the sys- 
tem size and the same behavior should be observed in 
the thermodynamic limit. We shall come back to this 
important issue in section [3 Another consequence of 
the definition of the field as an extremal quantity is that 
it exhibits large sample-to-sample fluctuations. In fact, 
it is found that the variance does not decrease with N, 
which means that each sample behave differently, even 
in the thermodynamic limit. It is however instructive to 
study in detail the probability distribution P^(H; to) for 
different values of to above and below a c . 

The evolution of the normalized histograms as a func- 
tion of system size is shown in Fig. One can see that 
the distributions are wide and rather complicated. On 
the one hand, there is a well-defined peak on the right- 
hand side of the histograms whose height increases and 
width decreases as N increases. This peak, however, does 
not exist for er < a c when to is in the range of the infinite 
avalanche (for instance m — 0.68 in the left panel of the 
figure) . On the other hand, there is another contribution 
which extends over a finite range and which is almost 
size- independent: it is responsible for the fact that the 
field is not self-averaging. By analyzing the sequence of 
microscopic states along each M-driven trajectory, we 
have checked that these two contributions come from the 
stable and unstable states, respectively. Since they are 
no other stable states than those belonging to the H- 
driven magnetization curve in the thermodynamic limit 
(as shown in Fig.[3Jl, we therefore conjecture that the dis- 
tribution PN(H;m) has the following asymptotic form: 
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FIG. 9: (Color on line) Normalized histograms of the output 
field H for selected values of m and different system sizes. 
The data correspond to 10 4 disorder realizations. 



Poo {H; to) = Q{m)5{H — H) + [l — Q{m)]w{H; to) (16) 

where S(H) is the Dirac function, H{m) is the field along 
the magnetization curve (i.e. the field taken as a func- 
tion of the magnetization), and w(H; m) is a continuous 
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distribution on a finite interval {H m i n (m), H max (m)}. 
Moreover, there is strong numerical evidence that 
H max (m) = H(m) for a > a c and for a < a c out- 
side the range of the macroscopic avalanche, whereas 
H max (m) < H(m) inside the infinite avalanche. 

The statistical behavior of the field for a given value 
of the magnetization is thus different above and below 
<j c . For a > er c , the most probable value of H(m) is the 
one corresponding to the ii-driven magnetization curve 
(the two protocols thus give the same field-magnetization 
diagram), but there is a finite probability that it takes 
a smaller value. For a < a c and m inside the range of 
the infinite avalanche, one has Q(m) — and the delta 
peak disappears. In this case, the value of the field is 
impredictable inside a finite interval. 



IV. RESULTS FOR THE CUBIC LATTICE 

Very similar results are obtained on the 3-D cubic 
lattice. In this case, however, the //-driven behavior 
cannnot be treated exactly in the thermodynamic limit 
and one must also perform simulations on finite systems. 

Fig.llOlshows the fraction of stable states along the M- 
driven trajectory. The curves are similar to the ones dis- 
played in Fig.0]except for a stronger asymmetry. Again, 
one finds that Q(m) — when m is in the range of the 
ii-driven macroscopic avalanche below <j c . 

Fig. ^2 shows that the internal energy obtained with 
the M-driven algorithm is still a self-averaging quantity 
both above and below a c (<r c ~ 2.2). However, it seems 
that the variance decreases slower than 1/L 3 when m 
is in the range of the infinite avalanche, a behavior also 
observed Ref.0] although the definition of the field is 
quite different. 

The comparison between the two algorithms for the 
average internal energy as a function of m is performed 
in Fig. El Note that there is again a double well struc- 
ture at low disorder but the two minima are very close to 
m = ±1 and hardly visible on the figure. There is only 
one minimum below a « 3, a value quite different from 
<j c . It seems again that the two algorithms give the same 
energy in the thermodynamic limit (outside the infinite 
avalanche) but the finite-size effects are more important 
than on the Bethe lattice. Finally, the histograms of the 
field H{m) are shown in Fig. 1131 The overall behavior is 
similar the one displayed in Fig.[§]and we thus conjecture 
that the asymptotic form of the probability distribution 
is given by Eq. (18). Note however that the continu- 
ous part w(H ; m) has a very different shape than on the 
Bethe lattice. 



V. DISCUSSION 

In the present paper, we have proposed a simple mod- 
ification of the standard single-spin-flip algorithm to 
study the magnetization-driven RFIM at T = 0. The 




FIG. 10: Fraction Q(m) of stable states along the M-driven 
trajectory on a cubic lattice (only the ascending branch is 
shown). The lattice size is L = 30 and the average has been 



taken over 3 x 
the eye. 



10 disorder realizations. Lines are guides to 



dynamics consists in flipping the spins one by one, choos- 
ing the spin with the largest local field. This allows to 
perform a detailed comparison with the microscopic tra- 
jectory of the system in the iJ-driven situation, above 
and below the critical disorder. It turns out that the 
two trajectories share the same metastable states in the 
thermodynamic limit, and we have computed the aver- 
age fraction of these states. An exact expression of this 
quantity has been obtained in the case of the Bethe lat- 
tice. Numerical simulations show that the two dynam- 
ics yield the same internal energy for a given value of 
the magnetization outside the range of the macroscopic 
avalanche. On the Bethe lattice, inside the macroscopic 
avalanche, the energy obtained with the M-driven algo- 
rithm also coincides with the one that can be computed 
analytically in the ii-driven case, using the solution of 
the self-consistent equations that describes the unstable 
branch of the hysteresis loop. 

The M-driven field-magnetization diagram exhibits 
some peculiar and annoying features that are due to our 
definition of the output field H as AU/AM: i) all closed 
loops have zero area, implying that there is no dissipation 
in the system; ii) H strongly fluctuates with m and these 
fluctuations are independent of the system size; iii) the 
sample-to-sample fluctuations of H also do not decrease 
with the system size. We have shown that this problem is 
related to the presence of a continuous part in the prob- 
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FIG. 12: Average internal energy per spin on the cubic lattice 
below and above o c . The symbols represent the results of the 
simulation of a system with size L — 30 with the //-driven 
and M-driven algorithms, as indicated. Averages are taken 
over 3 x 10 4 disorder realizations. Lines are guides to the eye. 



FIG. 11: (Color on line) Variance V u (m) of the internal en- 
ergy per spin u(m) on the cubic lattice for selected values of 
m as a function of system size. Averages are performed over 
3 x 10 4 disorder realizations. For the sake of clarity, the vari- 
ances for m = —0.75, m = 0.75 and m = 0.95 are divided by 
10, 100 and 1000 respectively. The lines are fits to the form 
V K (m)~L- p . 



ability distribution of H, which corresponds to the field 
associated to the unstable states. We now discuss some 
possible modifications in the definition of the dynamics 
or of the field. 

The first one is to allow for an additional relaxation 
of the system using the Kawasaki dynamics, which is the 
standard dynamics for a situation with a conserved or- 
der parameter. Specifically, one could imagine to first 
flip the spin with the largest local field (so to change M) 
and then perform all possible exchanges between nearest- 
neighbor spins of opposite sign that decrease the energy. 
This procedure is certainly more in the spirit of the lo- 
cal mean-field calculations that have been performed in 
Ref.Q at finite temperature. Even without changing the 
definition of the field, one may hope that the fluctuations 
of H with M will be weaker and will perhaps decrease 
with N. However, preliminary simulations show that the 
energy of some states cannot be decreased by exchanging 
nearest-neighbor spins and that many of the final states 
are still unstable with respect to the (Glauber) single- 
spin-flip dynamics. It would be interesting to perform an 



extensive study in order to understand if this behavior 
changes when increasing the system size. On the other 
hand, it must be emphasize that the stable states are now 
different from the ones visited by the //-driven dynamics. 
Moreover, this procedure does not solve the problem of 
the definition of the field (and the correponding absence 
of dissipation). 

A second possibility is to keep the same dynamics as 
in this work, but to change the definition of the output 
field. Indeed, it is very likely that the magnetization (or 
any other extensive variable) can only be controlled at 
the macroscopic level within a certain resolution. For in- 
stance, in an experiment performed at a constant rate 
dM/dt, one probably measures not the instantaneous 
force (in the present case — /™ aa: (m)) but some average H 
over a certain range Am (which could even depend on the 
driving rate). In this case, one can easily check that all 
fluctuations are suppressed in H in the thermodynamic 
limit since imposing a fixed resolution Am implies to 
take averages over larger and larger intervals AM when 
increasing N. H is then also self-averaging. One may 
also imagine that the apparatus that measures the field 
(or the force) cannot ajusts itself to the force infinitely 
fast or that there is some threshold value. Of course, in 
all these cases, the results are machine-dependent. Care- 
full "M-driven" experiments with different set-ups and 
different driving rates are thus needed in order to better 
resolve these issues. 

Finally, from a theoretical point of view, one cannot 
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FIG. 13: (Color on line) Normalized histograms of the output 
field H for selected values of m and different system sizes. The 
data correspond to 3 x 10 4 disorder realizations of a system 
with size L — 30. 



discard the possibility that there does not exist any sat- 
isfatory definition of the output field when using Ising 
variables. An alternative approach using continuous vari- 
ables has been proposed in Ref.0. 
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